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Energies and wave functions are calculated for d-wave quasiparticles in the mixed state using the 
formalism of Franz and Tesanovic for the low-lying energy levels. The accuracy of the plane-wave 
expansion is explored by comparing approximate to exact results for a simplified one-dimensional 
problem, and the convergence of the plane wave expansion to the two-dimensional case is studied. 
The results are used to calculate the low-energy tunneling density of states and the low-temperature 
specific heat, and these theoretical results are compared to semiclassical treatments and to the 
available data. Implications for the muon spin resonance measurements of vortex core size are also 
discussed. 
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I. INTRODUCTION 



The nature of the low-lying excitations in the mixed 
state of a <i-wave superconductor is both an interesting 
quantum mechanics problem and important for under- 
standing the behavior of high temperature superconduc- 
tors in a magnetic field |5]-k|. Volovik jl| first studied 
this problem in the semiclassical limit, where the d-wave 
quasiparticles are Doppler shifted by the local superfluid 
density. The shifting of quasiparticle energies results in a 
non-zero density of states at zero energy proportional to 
the square root of the magnetic field. Volovik's solution 
has been applied to calculations of the specific heat 
thermal conductivity PJT(i|], and nuclear magnetic relax- 
ation rates jLl 12|. It has motivated useful discussions 
of the scaling behavior of the specific heat by Simon and 
Lee JDj, Kopnin and Volovik ]14| ], and others. 

The presence of a magnetic field and its associated vor- 
tex lattice affects the motion of quasiparticles in four 
distinct ways. First, the quasiparticles, which carry cur- 
rent, move in the magnetic field which is approximately 
uniform for an extreme type-II superconductor. Second, 
although the field is approximately uniform, it is not ex- 
actly so, and therefore the quasiparticles experience mag- 
netic field gradients. However the direct effect of these 
gradients is rather small. Third, there are supercurrents 
associated with the curl of the field, and the quasiparti- 
cle energies are affected by the corresponding superfluid 
velocity through which they move. For a uniform super- 
fluid velocity field, the effect would be a simple Doppler 
shift of the energies. However, for inhomogeneous super- 
fluid velocities the problem is more complicated. Fourth 
and finally, the magnitude of the superconducting order 
parameter is inhomogeneous in the mixed state, although 
this is mainly apparent within a coherence length of each 
vortex core where this magnitude falls to zero. For an 
extreme type-II superconductor, this represents a very 
small fraction of the sample for fields well below H C 2- 

Volovik's approach neglects the magnetic field and its 
gradients as well as the inhomogeneous order parameter 
amplitude and focuses only on the local supercurrent ve- 
locity. It assumes that the quasiparticle wave function 



can be thought of as a wave packet which is localized in 
a region over which the magnitude and direction of the 
superfluid velocity are relatively uniform. The energy 
of a low- lying, d-wave quasiparticle depends linearly on 
q = k — k a , where k a is the wave vector of the nearest 
node. If a quasiparticle is localized in a region of size 
smaller than l/\q\, then the spread in its energy will be 
larger than its average energy and the wavepacket pic- 
ture does not work. For the superfluid velocity to be 
relatively uniform in a region, the size of the region must 
be smaller than the distance to the nearest vortex core 
and certainly smaller than, d, the distance between vor- 
tices. Let us apply the above considerations to the lowest 
energy quasiparticle band. This corresponds to a quasi- 
particle, near node u, with wave vector q perpendicular 
to k v localized in a region of size £. For the wavepacket 
picture to apply, the energy, E, of the quasiparticle must 
be greater than Uv^k jl, where «a is the quasiparticle ve- 
locity along q. However for the superfluid velocity to be 
uniform in the region of size £ it is necessary that £ <C d. 
Combining these two conditions we obtain the require- 
ment that E ^> Tiv&ir/d. For energies less than this the 
wavepacket picture breaks down and a full quantum me- 
chanical picture is needed. This energy range is readily 
accessible via specific heat measurements below 10 K in 
fields of one to several Tesla. It is this energy region that 
is the main focus of this paper. 

Recently, Franz and Tesanovic Ji| (FT) have derived 
a quantum mechanical theory of the mixed state of a 
d-wave superconductor, which involves a singular gauge 
transformation that maps the original problem of su- 
perconducting quasiparticles in a magnetic field onto an 
equivalent one of quasiparticles in a periodic potential. 
The latter problem may be solved using conventional 
band-structure methods. 

In this paper, we investigate the low-energy proper- 
ties of a d-wave superconductor in the mixed state us- 
ing the theory derived by Franz and Tesanovic The 
most direct experimental probes of these properties are 
the low-energy tunneling density of states and the low 
temperature specific heat. In order to calculate these 
quantities reliably, we have investigated the numerical 
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problem of Dirac quasiparticles in the periodic poten- 
tial of the vortex lattice, focusing on the simplifications 
that result from the fact that the anisotropy of the Dirac 
cones, <xp> = vp/v& ^> 1. As discussed by Mel'nikov 
Jl5|] , such large anisotropy makes the problem approxi- 
mately one-dimensional. Mel'nikov described how to ob- 
tain solutions to the one-dimensional problem, but he 
then confined his analysis to the semiclassical versions of 
these solutions. We have explicitly evaluated the quan- 
tum mechanical solutions in this one-dimensional limit 
and used them as a test of the accuracy of approximate 
plane-wave solutions. We then show how to improve on 
the one-dimensional solutions by including a small num- 
ber of plane wave basis states for the transverse direction, 
and we study the convergence of this approach. 

The remainder of this paper is organized as fol- 
lows. Section II addresses the computational problem 
of calculating quasiparticle energies in the lowest bands 
in the magnetic Brillouin zone, comparing exact and 
plane-wave-expansion solutions for the simplified one- 
dimensional problem and then comparing 1-D and 2- 
D plane-wave expansion solutions for various choices of 
plane-waves. Section III presents results for the local 
tunneling density of states, Sec. IV re-interprets recent 
muon spin resonance measurements of the vortex core 
size in terms of a scaling picture, and Sec. V presents 
calculations of the low temperature specific heat, com- 
paring them to predictions based on Volovik's approach 
and to experimental data. 



II. THE COMPUTATIONAL PROBLEM: 
CALCULATING THE ENERGIES IN THE 
LOWEST BANDS IN THE MAGNETIC 
BRILLOUIN ZONE 



The quasiparticle wavefunctions are described by the 
BdG equations, Hip = eip where ip = (u(r), v(r)) T , and 



n = 



n e a 
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with H e = (p — - A) 2 — ep. The gauge invariant form of 
the gap operator, A, for a <i-wave superconductor can be 
written as (see Ref. p6i for the details) 
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where, for notational convenience, we have chosen to 
orient our axes along the directions of the gap nodes, 
at an angle of 7r/4 with respect to the orientation of 
the Cu02 planes, pf is the Fermi momentum, and 
A(r) = | A(r)\e l(p ^ is the Ginzburg-Landau order pa- 
rameter. Since we are working in the intermediate field 
regime (H c \ < B « H C 2) of an extreme type II super- 
conductor, we can assume that the magnitude of the gap 
is constant everywhere, except at the vortex cores, and 
that the magnetic field distribution and local superfluid 
velocity can be described by the London model Jl7|] . 



In order to diagonalize the Hamiltonian in Eq. ([j]) one 
would like to remove the order parameter phase from the 
off-diagonal components of H. It is desirable to choose a 
transformation which is both single- valued and treats the 
particles and holes on an equal footing. This is accom- 
plished by the bipartite, singular gauge transformation 
of FT: 1 



H -> U~ l HU, U 
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where <PA( r ) + ^s(r) = <f(r), and ^^(r) and <y9s(r) are 
the contributions to the order parameter phase from the 
A and B sublattices of the vortex lattice. The sublattices 
are chosen so that there are equal numbers of A and B 
vortices, with two vortices per magnetic unit cell of the 
vortex lattice. The vortex lattice configuration analyzed 
in this paper is shown in Fig. [l| Note that the fact that 
the x- and y-axes of the A and B sublattice unit cells 
are oriented along nodal directions means that nearest- 
neighbor lines of vortices are oriented along the a- and 
6-axes of the underlying crystal lattice. 
Under this transformation H becomes 
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where 
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with the superfluid velocities 

<(r) = - (hV^ - - A) , n = A, B. (6) 

m \ c / 

Note that v^(r) + vf (r) = 2v s (r). Since the vortex lat- 
tice is periodic, the superfluid velocities can be written 
as Fourier sums 
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where K — (2tt /d)(m Xl m y ), d = ^j2<S>a/B is the size of 

the magnetic unit cell, and 5^ = ±(<i/4, d/4) is the dis- 
placement of A or B vortices from the center of the unit 
cell (see Fig. |l|). 

Linearizing the Hamiltonian in Eq. (^) at the node 
k = (kp,0) we find that Hab — Ho + H' with 
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and 
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Aq/pf is the 



where vf is the Fermi velocity, and va 
slope of the gap at the node. 

At the node k = (kp,0) the free Dirac Hamiltonian 7i 
has the familiar Dirac cone spectrum 



^k 2 x + k 2 y /a 
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where {k x ,k y ) = at the node. The quasiparticle mo- 
mentum along the nodal direction is k x ~ e/hvp with a 
corresponding wavelength of X x ~ Tivf/s- If the energy is 
low enough (e < hvp/r x ), k x will be confined to the first 
magnetic Brillouin zone (see Fig. ^ and the wavelength 
X x will exceed the intervortex distance r x , crossing the 
boundaries of several unit cells of the vortex lattice. For 
large values of the anisotropy ajj the momentum parallel 
to the Fermi surface k y ~ aok x will be much larger than 
that along the nodal direction. The quasiparticle wave- 
function will thus be localized in the direction parallel to 
the Fermi surface, but will be extended and will feel the 
average effect of the superfluid velocity fields of several 
vortices along the nodal direction. 

Since the potential TL' is periodic we can expand the 
quasiparticle wavefunction in the plane-wave basis: 



*k(r) =J2 e 



i(k+K)-r 



V^ K (k). 



(11) 
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The periodic potential of the vortex lattice will be respon- 
sible for the interaction of the Fourier components ^K(k) 
and V'K' (k) . If we are interested only in energies below 
a cutoff energy, E c , for which the momentum k x lies well 
within the first MBZ we can make the approximation 
that the quasiparticle wavefunction is one-dimensional 
and ignore the interaction of those Fourier components 
that are at different values of K x . We therefore write 



*k(r) 



■ e i{ k v+ K v)y 



<K(k). 
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If, however, E c is high enough that k x exceeds the bound- 
aries of the first MBZ (see Fig. ||) we can make the 
assumption-since the Fourier sum is dominated by com- 
ponents whose values of K are bounded by the constant 
energy contour at _E c -that 



*k(r) 



K a K c /a D 



V>K(k), 
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where K c is the cutoff wave vector along the y direction. 

Such plane-wave expansions can be computed numer- 
ically to obtain the excitation spectrum for the quasi- 
particles in a periodic vortex lattice. The solution to 
the problem using Eq. ( |TT| ) has been studied in detail by 
FT j|, whereas Marinelli, Halperin and Simon (MHS) 
[jl8| studied solutions to Hab defined in Eqs. (||) and 



(g) in position space. Both groups found that the con- 
vergence of the plane-wave expansions was slow. Since 
we are specifically interested in the low energy and low 
temperature properties which are largely determined by 
the lowest band of the excitation spectrum, we will focus 
next on obtaining an analytical solution to the linearized 
Hamiltonian with the approximation that, for large a_o, 
the quasiparticle wavefunctions are approximately one- 
dimensional. Having obtained both analytical and nu- 
merical solutions to this one-dimensional problem, we 
will then examine how adding more transverse wave vec- 
tors , as in Eq. (|l3|), allows us to approach the exact 
numerical two-dimensional result, using Eq. (pT|). 



A. The 1-D analytical solution 

At low energies and for large values of old the wave- 
functions are localized in the y direction and extended 
along the x direction. jl5| This suggests the following 
basis as a useful starting point: 



(14) 



A",, 



As we shall see, the Fourier components ip(k x + K x ,y), 
for different K x , are spatially well separated in the y di- 
rection. Their interaction is consequently negligible, and 
we can assume that the Hamiltonian is diagonal in the 
quantum number K x . This allows us to replace the pe- 
riodic potential Tt' in Eq. (j^j, which in principle scatters 
quasiparticles between states with different values of K Xl 
with its effective potential averaged in the x direction 
which is diagonal in K x . The result is 
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where Titp — eip, e = E/Eq where Eq = 2Trhi)F/d, 
z = y/d, and $(z) = z — (n + i) where n < z < n + 1. 
Note that in these units n — \ < q x < n + 5, where 
q x = (d/2Ti)(k x + K x ). The potential along the y direc- 
tion consists of two periodic sawtooth functions with dis- 
continuities lying along the averaged vortex lines of the A 
and B sublattices. At sufficiently low energies the quasi- 
particles will be bound in the y direction by the poten- 
tial barriers which lie at the discontinuities in $(z). Our 
picture is thus one of quasiparticles that travel as plane- 
waves along the nodal direction but are bound within 
potential wells-created by the averaged vortex lattice-in 
the direction parallel to the Fermi surface. Note that, 
at low energies, the Fourier components tp(k x ,y) and 
ijj(k x + K x , y) have negligible overlap, as they lie in sep- 
arate potential wells along the y direction. 
By making the substitution 
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Eq. [l5| can be rewritten as 
a dz 



e + a x [ g£ + 



<p[z) = 
(16) 



where a = Itxcxd and &i are the Pauli matrices. The 
function 

= *(*--) +*(* + -) 

is a sawtooth with a slope of +2 and a period of 1/2, and 
the function 

* 2 {z) = 9>{z-\)-$(z+±) 

is a step function which oscillates between —1 and +1 
with a period of 1 . Since the potential is periodic we can 
solve Eq. (|l6|) within a unit cell and use Bloch's theorem 
to extend the solution over all of z. The solution within 
a unit cell (see Appendix [a| for the details) is given in 
terms of the parabolic cylinder functions D p (z): |l9| 



(17) 



for n — j < z < n + j and 
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with zi = \fa[z — e — n) and 

^2 = \/a(z — e — n — 5 ). 

These solutions can be matched at the boundaries of 
the unit cell (see Appendix B) to obtain an exact excita- 
tion spectrum for the one-dimensional, averaged Hamil- 
tonian. The resulting spectrum for anisotropy old = 7 is 
shown in Fig. ||. It is useful to note that the energy scale 
Eq is approximately given by E sa 185Vi3 T -1 / 2 K. 



B. Comparison of 1-D analytical and plane- wave 
expansion results 

Using the plane- wave expansion of Eq. (Jl^) we can nu- 
merically diagonalize the Hamiltonian to obtain an exci- 
tation spectrum that can be compared with the analytical 
results. To the numerical accuracy of the diagonalization, 



these two methods yield identical results for the disper- 
sion along k x as shown in Fig. |§|, where 61 reciprocal 
lattice vectors (RLVs) have been kept in the plane-wave 
expansion. |2(it| The dispersion along k y calculated from 
the 1-D plane-wave ex pan sion is also shown in Fig. [jj. 
As discussed by MHS, [|l8| the dispersion away from the 
r point along k y is more strongly renormalized by the 
supercurrents, leading to an enhanced effective old- For 
old > 10 there is essentially no discernable dispersion 
along k y for the lowest bands (as calculated in either the 
1-D or full 2-D plane- wave expansions), further suggest- 
ing the validity of a one-dimensional approximation. 

Since both the energy and momentum axes scale as 
\f~B these spectra apply to all values of the magnetic field 
within H c i <<£?<< H C 2- As the anisotropy increases, 
the gap between the lowest band of the spectrum and 
the E = axis quickly narrows. At ao = 14 (Fig. |^) the 
spectrum is close to forming a line quasinode, in agree- 
ment with the results of FT g, and of MHS |lf], which 
suggest that a line quasinode first appears at ajj ~ 15. 



C. Comparison of 1-D and approximate 2-D 
plane-wave calculations 

The results of numerical diagonalization calculations 
of the excitation spectrum of the quasiparticles in the 

1- D averaged potential and the exact 2-D potential at 
different values of the anisotropy are shown in Figs. |] 
and ||. The 1-D spectra show good qualitative agreement 
with the 2-D spectra, capturing the major features of the 
lowest bands, including the line quasinodes that form at 
large values oi an- However, the 1-D treatment is un- 
able to accurately represent, quantitatively, the behavior 
of the full 2-D spectrum. In particular, as can be seen 
from Figs. [| and the 1-D approximation cannot be 
used to quantitatively determine the size of the minigaps 
which lie along the line quasinodes. An analysis of the 
TD spectra for several values of au shows that the size 
of the smallest minigap at fixed old is S g cx e ~ maD where 
to w 0.18. Unfortunately, the slow convergence of the 

2- D reciprocal lattice sums-due to the divergence of the 
superfluid velocity at the vortex cores (discussed in more 
detail by Vafek et al. jl6|])-makes it very difficult to ac- 
curately determine the size of these minigaps in the full 
2-D calculation. 

Nonetheless, we believe that the 1-D treatment (which 
is far less computationally intensive than the 2-D prob- 
lem) captures and elucidates the important physics of 
the lowest bands of the quasiparticle excitation spectrum 
and is therefore a useful tool that helps us understand the 
physical behavior of the quasiparticles in the mixed state. 
In particular we will use the 1-D energies and wave func- 
tions to calculate the local tunneling density of states and 
the specific heat. 

Next we compare the results of the 1-D calculations to 
finite 2-D plane wave expansions using a grid of N x x N y 
reciprocal lattice vectors. For example, in Fig. || re- 
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suits are shown for ajj = 14, comparing the 1-D case, 



1, N y = 41, to N x = 5, 9, 13, 21, and 29, N y = 41. 



Similar results are shown in Fig. ^ for — 20 and 
= 61. One of the most striking features of both fig- 
ures is the complete inscnsitivity of the linear branch near 
the r point to the number of plane-waves in the calcu- 
lation. For this branch, it appears that the 1-D energies 
are essentially exact. For other low-energy branches and 
general points in the Brillouin zone, the plane-wave ex- 
pansion seems to converge smoothly. The only patholog- 
ical behavior occurs near the quasinodes, where both the 
positions and the values of the minima converge slowly. 

III. LOCAL TUNNELING DENSITY OF STATES 

In this section we show results of calculations of the 
local tunneling density of states (TDOS) of the quasipar- 
ticles in the lowest band of the energy spectrum using the 
one-dimensional plane- wave expansion of Eq. (|l2|). The 
TDOS is: |^|| 

N(r, 15) = - A £ K^rjf/W - E) 

+ |^k,^(r)| 2 /'(£k, M + S), (19) 

where f'{x) is the derivative of the Fermi function, k is 
the set of wave vectors in the magnetic Brillouin zone, /j, 
is the set of energy bands (restricted in this case to the 
lowest positive and negative energy bands) and v is the 
set of four Dirac nodes. The normalization factor 
is equal to the number of spins divided by the number 
of wave vectors in the magnetic Brillouin zone. The ex- 
tra factor of 2 is an artifact of the normalization in the 
diagonalization routine used. 

The plane-wave expansion was done at the node k = 
(fejr, 0_V It is easy to show that by taking e^,^ — > — £k )( u hi 
Eq. (M) one obtains the contribution from the opposite 
node at k = (— fcp,0). Within the 1-D approximation, 
these two nodes give the y dependence of the TDOS, 
and the other two nodes at k — (0, ±fcp) give the x de- 
pendence. The y (or x) dependence of the TDOS at 1 
Kelvin and at a field of 1 Tcsla, for two different values 
of the anisotropy ao, and at three different energies is 
shown in Figs. ^ and ^. 

One can see that the TDOS has the periodicity of, and 
is sharply peaked at, the vortex lines. The TDOS falls 
to a broad minimum in the regions between the vortices. 
The shoulders on either side of the peaks come from the 
states within the lines of quasinodes that form at large 
values of ao- At ao — 20, the size of the gap in the line 
quasinode has decreased and a second line quasinode has 
started to appear (see Fig. ||). Both these features con- 
tribute to the very distinct shoulders on either side of the 
peak in the TDOS in Fig. |. 

Figure [lO] shows the zero-bias two-dimensional TDOS 
as a sum over the four nodes. This result is in qualita- 
tive agreement with the semiclassical calculation of the 



TDOS by Mel'nikov |2J. The vortex lattice geometry of 
our paper is, in Mel'nikov's notation, a Type II lattice 
with a — 1/2. This gives a TDOS that is proportional to 



d/2 



(20) 



where $(z) = 2z — (2m + 1). The semiclassical TDOS 
of Mel'nikov thus has the profile of a triangle wave along 
the x and y directions. The fully quantum mechanical 
results shown here follow this profile, but exhibit addi- 
tional structure that arises from the quasiparticle states 
near the quasinodes. 

We note that only half of the bright spots in Fig. [h] 
lies at vortex positions while the other half lies halfway 
between vortices. For example, in Fig. [To] the bright 
spots at the corners and at the center of the figure might 
correspond to vortex sites. The other bright spots are 
then the result of the overlap of the sharply peaked tun- 
neling density of states which extends from each vortex, 
parallel to the four node directions. It is an artifact of 
the 1-D model that, for the case of a square lattice, these 
overlaps have a peak tunneling density of states equal to 
that of a vortex core. This artifact is less evident in more 
general, centered rectangular lattices or, in particular, for 
the hexagonal lattice. p5|,p3| 



IV. MUON SPIN RESONANCE 



Two important simplifying assumptions in this model 
are that the superconducting coherence length is negligi- 
ble and that the penetration depth is large compared to 
the distance between vortices. As a consequence of these 
assumptions, the intervortex spacing is the only length 
scale in the problem. This allows us to present results 
scaled to this length as is done above for the tunneling 
density of states. 

In addition to the tunneling density of states, one could 
also use the wave functions generated by these calcula- 
tions to compute the pattern of the two-dimensional su- 
percurrent density. This would, of course, not be a self- 
consistent result, but it would be an improvement over 
the initial form for the supercurrent density correspond- 
ing to Eq. (7). Without actually doing this calculation, 
we know that the resulting pattern would be a function 
of r/d and hence that all lengths would scale as 1/y/B. 

This picture, in which the vortex lattice constant pro- 
vides the only length scale, is supported by the self- 
consistent calculations of Franz and Tesanovic for a single 
d-wave vortex. |24| In Fig. 1 of Ref. |24|] and the accom- 
panying discussion, it is shown that, for systems with 
very short coherence lengths, the spatial dependence of 
the gap function outside the core has a scale-independent 
power-law dependence, approaching its asymptotic value 



roughly as 1/r 
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The above discussion provides a natural explanation 
of the muon spin resonance results of Sonier and co- 
workers [ p5| who found that the vortex core radius, de- 
fined as the radius at which the supercurrent density 
has its maximum, grows large at low field. In fact, an 
excellent fit to their data can be obtained by assum- 
ing that the vortex core radius scales as \j\[B, as is 
shown in Fig. |H| The coefficient of 1/yB from the fit is 
r = 46.3 ± l.hAT 1 / 2 j\fB. Since the vortex lattice con- 
stant d for the A or B sublattices is d = 632AT 1 ^ 2 /y f B, 
this maximum occurs at about 7% of d or equivalcntly 
at about 10% of the intervortex spacing. It would be in- 
teresting to test this result at higher fields to see if this 
scaling breaks down and if r$ saturates at a constant 
value limited by the coherence length £o as one might 
expect. 



where No is the zero-field density of states. For nonzero 
E we find that 



N(e) = N(0) 



6eVl - 4e 2 + (8e 2 + l)sin- 1 (2 £ ) 



8e 



for < |e| < 1/2, where e = E£/hv F , and 



A^) = fiV(0)( £+ i 



, (25) 



(26) 



for |e| > 1/2. Note that this is the contribution to the 
total density of states for 2 spin states from one of the 
four nodes. 

A more realistic calculation of the semiclassical DOS 
can be made for a square vortex lattice if we write the 
superfluid velocity as the Fourier sum 



V. THE DENSITY OF STATES AND THE 
SPECIFIC HEAT 



A. Semiclassical DOS 

We start by calculating the density of states for the 
semiclassical (SC) approximation, in which the energy is 
Dopplcr shifted by the local superfluid velocity v s (r) : 



E(k, r) = hk F x • v s (r) ± ^(hv F k x ) 2 + (hv A k y ) 2 , (21) 

where the spectrum has been linearized around the node 
k = (fc F ,0). 

The local superfluid velocity far from the vortex is 
v s (r) = (h/2mr)4>. In the commonly employed "single- 
vortex approximation," the associated density of states 
is 
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where Q = 2Tt(mx + ny)/a and a = ^/&q/B. Note that 
we are now orienting the x and y axes along the nearest- 
neighbor directions of the square vortex lattice. The cor- 
responding density of states, 
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can then be calculated numerically using this more ac- 
curate expression for v s . The semiclassical density of 
states, as calculated for both the single-vortex approxi- 
mation and the square vortex lattice, is shown in Figs, [l^ 
and [l3| One can see that the square-lattice DOS is about 
30% lower at zero energy than that calculated using the 
single- vortex approximation. This lowering is caused by 
the disappearance of v s (r) at high symmetry points on 
the vortex-lattice unit-cell boundary. [p6| 



where the factor of 2 accounts for spin degeneracy. V/ w 
is the total area of the CuO planes in the sample, where 
V is the volume of the sample and w is the average sep- 
aration between the planes, and ir£ 2 = &o/B is the area 
of one unit cell of the vortex lattice. The integral is over 

€ = ^{hv F k x ) 2 + (flVAky) 2 - 

In the absence of a magnetic field, with no Doppler 
shift, 



No(E) 



V 



irh v f vaw 



\E\ 



(23) 



Putting the magnetic field back in, the density of states 
has the intercept 



JV(0) = -N 



hv F 



(24) 



B. Quantum DOS 

The quantum density of states is calculated from the 
quasiparticle energy spectrum at the node k = (k F ,0): 



N(E) 



V 



d 2 N k w 



Y,S(E~E nk ), 



(29) 



where 
tor in 
counts 
sional 
byl = 
tors in 
planes 



n labels the energy bands and k is a wave vec- 
the magnetic Brillouin zone. The factor of 2 ac- 
for spin degeneracy. In order to clarify the dimen- 
analysis, we have multiplied the usual expression 
V / (voN^d 2 ), where is the number of wave vec- 
the MBZ, and V/w is the total area of the CuO 
in the sample. The energy in this expression is in 
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units of 2irhvF/d. In order to compare this result with 
the semiclassical result we simply write N(E) in units 
of N Q (hv F /£) [see Eq. ||l noting that I = d/V2n. Re- 
sults are shown in FigsTp.2j and |l3j, where comparison is 
made to both the SC single-vortex approximation and 
to the SC square-lattice DOS. Note that both axes scale 
as l/£ oc v~B. The dotted line shows the commonly em- 
ployed single vortex SC DOS to be roughly twice as large 
as the quantum 1-D DOS in the low-energy region. The 
quantum DOS rises more quickly with energy and the SC 
and quantum DOS match up at higher energy andare in- 
distinguishable for energies above 3Ey . The discrepancy 
between the SC square-lattice calculation and the quan- 
tum DOS at low energies is due to quantum effects that 
average over the rapid variations in the direction of v s (r) 
near the vortex cores as well as near the high symmetry 
points on the unit-cell boundary. Of course, disorder ef- 
fects on the vortex lattice and the quasiparticle energies 
will also affect the average magnitude of the low-energy 
DOS in both the SC and quantum cases. [|]j2(| 

The 1-D calculation of N(E) for a D = 20 (Fig. |TJ) is 
in good agreement with the corresponding 2-D calcula- 
tion of FT [||], reproducing all of the major features at 
low energies. The overall magnitude of the 1-D DOS is 
slightly reduced from the full 2-D calculation. The 1-D 
calculation, by essentially averaging in one direction, un- 
derestimates the effect of the supercurrents, which push 
states to lower energies, as can be seen from the band- 
structures shown in Figs. ^ and |^. The full 2-D quantum 
DOS in Ref. j§, is about 10% higher in magnitude than 
the 1-D approximation, but is still noticeably lower in 
magnitude than the SC square-vortex-lattice result. 



C. Scaled C V (T,B) 



The heat capacity of a fermion gas is 



C = 2l3k B J2 



df(E k ) 
dE k 



El 







N T {u/[3) 



1 + cosh 



-du . 



This is the expression used to calculate the specific heat 
(C v = C/V) from the total density of states N T (E). 
The total density of states in Eq. ( |30| ) is a sum over the 
density of states for one spin at each of the four nodes. 
ThusJV T (^)= 2N(E) where N(E) is the semiclassical 
[Eqs. |2^ & ^6) or quantum mechanical [Eq. ^9| density of 
states calculated in the previous section. Therefore, the 
specific heat at constant volume is 



P T 

r - 1 B 



N(u/f3) 



1 + cosh u 



-du. 



(31) 



The specific heat for the 1-D calculation is shown for var- 
ious values of an in Fig. Again, both the C/T axis 



and the T axis scale as yB, in agreement with the gen- 
eral scaling predictions of Volovik [[[) , and Simon and Lee 
The C v /T is linear at higher temperatures, flattens 
out as the temperature is decreased, and then increases 
to a peak at even lower T before rapidly falling, with 
a tiny shoulder on the way down (see inset of Fig. fbi] ), 
to zero at T = 0. The large peak both sharpens and 
moves closer to the T — axis as the anisotropy ajj is 
increased. The behavior of this peak suggests that its 
presence is due to the low energy peaks in the DOS, par- 
ticularly the van Hove singularities that occur just above 
Ey = 0.25 for ac = 14 as these contribute significant 
weight to the DOS. A narrow peak at E in the DOS will 
typically show up as a peak in the specific heat near E/2. 
P9j Comparing the 1-D and 2-D dispersions and DOS, 
we expect this peak to shift to slightly lower energy and 
to sharpen in the full 2-D calculation of the specific heat. 

The SC specific heat, for the square lattice and for the 
single-vortex approximation, is shown in Fig. Es), along 
with the 1-D specific heat. The temperature is in units 
of E v = hvp /£ and C v is in units of 



(32) 




Again, the main difference between the SC and quan- 
tum specific heat is that the SC specific heat larger in 
magnitude. Both exhibit the same scaling with magnetic 
field and with an- The quantum specific heat exhibits 
structure at the lowest temperatures which is a reflection 
of the structure in the low-energy DOS. 

In order to make comparisons with experimental re- 
sults, we use the numbers in Chiao et al. |27| for YBCO: 
vf — 2.5 x 10 7 cm/s, ao — 14, and w = 5.85 A. The mo- 
lar volume of YBCO is V M = 104.38 cm 3 /mole. @ With 
these numbers we obtain an intercept for the SC single 
vortex calculation of 0.91 (mJmoP 1 )K _2 T -1 ^ 2 in seem- 
ingly excellent agreement with the experimental \f~B~ co- 
efficient of 0.91 (mJmor^KT 2 of Moler et al. §. How- 
ever, since this approximation overestimates the zero en- 
ergy specific heat by roughly a factor of 2, this agreement 
is fortuitous. The quantum specific heat for ac = 14 

(30) 

flattens out at approximately 0.5 (mJmol : )K 2 T ' . 



The Geneva group of Junod and co-workers has re- 
ported a number of results Jt|]3C|] for the specific heat 
of very high quality YBCO crystals, grown in BaZrOa 
crucibles and doped to O7.00 so as to minimize the ef- 
fects of oxygen chain vacancies. In the earlier of these, 
Revaz et al., J7j the vortex contribution to the specific 
heat was obtained by subtracting C(B _L c, T) from 
C(B || c,T), the idea being that both lattice and mag- 
netic impurity effects would cancel out in this subtraction 
and that the vortex contribution to the spceific heat for 
B _L c is small. In the more recent preprint by Wang 
et al. [|0), results are presented for C(B,T) - C(0,T) 
for B [|c. For our purposes, these data are more di- 
rectly useful since they involve only the single field di- 
rection that we have studied. Furthermore, the results 
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should be reliable, since the samples show very little sign 
of point magnetic impurities. Wang et al. |30t] find a 
T ^ intercept for (C(B, T) — (7(0, T))/T of 1.34 ±0.04 
(mJ mol-^K^T- 1 / 2 . 

In order to compare our theoretical results to the T 
and B dependence found by Wang et al., |}0) we need 
to subtract the specific heat in zero field from that in a 
field. The result is shown in the inset of Fig. [ll]. It is 
interesting that the structure that we find at the lowest 
temperatures could easily be attributed to Schottky-type 
anomalies in the data. In fact Wang et al. j30| show fig- 
ures with and without subtraction of an assumed Schot- 
tky anomoly, and the latter better resembles our theo- 
retical results. It is tempting to suggest that the experi- 
mentally observed low temperature structure in samples 
with the least magnetic impurities is actually due to the 
structure in the quasiparticle density of states. However, 
since the magnitude of the observed field-dependent spe- 
cific heat is more than twice as large as the calculated 
value, such detailed comparisons between theory and ex- 
periment are probably premature. The effect of disorder 
on the quasiparticles would likely increase the low energy 
specific heat, since it increases the low energy DOS. On 
the other hand, disorder in the vortex lattice may equally 
well decrease the low-energy specific heat by reducing the 
local supercurrent velocity. Therefore, it is not obvious 
that disorder, in itself, can account for the differences 
between theory and experiment. 

VI. CONCLUSIONS 

By making the approximation, in the mixed state, that 
the low-energy quasiparticle states in the Dirac nodes are 
essentially one-dimensional, we have been able to obtain 
analytical results for the quasiparticle wave functions and 
energy spectra. The 1-D approximation to the FT Hamil- 
tonian [|j elucidates the physics of the interaction of the 
quasiparticles in the lowest-energy bands with the vortex 
lattice: the quasiparticles travel as plane- waves along the 
directions of the gap nodes, and are confined by the pe- 
riodic potential of the vortex lattice in the direction of 
the Fermi surface. Using these exact analytical results, 
we were able to show that the approximate plane-wave 
solution for the same problem converges rapidly. The 
1-D approximation is able to reproduce the important 
features of the 2-D plane-wave expansion in the lowest 
bands. 

We have presented calculations of the tunneling den- 
sity of states, which are in qualitative agreement with 
the semiclassical results of Mel'nikov [^3| but which also 
show spatial structure due to the energy dispersion of the 
low-lying states. The density of states at zero energy for 
the quantum problem is significantly lower (by a factor of 
2) than the commonly-employed semiclassical result for 
a single circular vortex, although this simple approxima- 
tion overestimates the semiclassical density of states for 
a square vortex lattice configuration by roughly 30% at 



zero energy. Thus this reduction arises from two sources: 
the larger area of low superfluid velocity in the Abrikosov 
lattice, compared to the case of a single vortex in a cir- 
cular unit cell of radius £, and quantum averaging of the 
superfluid velocity for quasiparticles in the first magnetic 
Brillouin zone. 

The specific heat has been calculated from the DOS 
of the 1-D plane-wave expansion and found to exhibit 
structure at low temperatures that is not present in the 
semiclassical approximation. In addition, the magnitude 
of the low- temperature specific heat is reduced by quan- 
tum effects. Since the values of the specific heat measured 
experimentally, j^J^,^] for parameters vf and va taken 
from other experiments, are already larger than the semi- 
classical results, the disagreement in magnitude with the 
quantum results are even larger. One possibility is that 
the discrepancy is due to the effects of disorder which, 
to date, have not been included in any quantum treat- 
ment of the specific heat. Other possibilities are that the 
parameters are such that ao is substantially larger than 
is currently believed, or that there are additional low- 
energy states not accounted for in the disordered d-wave 
model but which exhibit similar magnetic-field depen- 
dence. 
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APPENDIX A: ANALYTICAL SOLUTION OF 
THE 1-D PROBLEM 

With the approximation that the potential is one- 
dimensional the quasiparticle Hamiltonian is 

(Al) 

The Hamiltonian can be rewritten as 

(A2) 

where $i(z) = $(2 - |) + <5>(z + j), and $2(2) = 
$(z — i) — <5>(z + |). Borrowing a trick from Mel'nikov 
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p5[ , one can insert I = 1/2(0";,; + a z )(&x + ®z) between 
7iiD and ifj an d then multiply Eq. ( |A2| ) on the left by 
l/2(o' a; + a z ). This transformation takes a x — > cr z and 
vice versa. We then write 



so that 



a dz 



ip(z) = l/2(a x +a z )i)(z), 



* 2 (z) 



/"(z) + a 2 ((z - c,) 2 - 6 2 + /(z) = 0, (A6) 

where j = 1, 2 and 6i = g™ + |, &2 = — j, Ci = e and 
C2 = e + 5 . If we let r,- = v^O 2 — c j ) then 



f"(r j )+(^-ab 2 j +i)f(r j ) = 0. 



(A7) 



^(z) 

= 0, (A3) 



This equation is solved by the parabolic cylinder func- 
tions (see Gradshteyn and Ryzhik Jigj] ) 



ffa)=D iXj [±(l+i)Tj], 



(A8) 



where a = 2naD- Writing ip(z) = {f{z) 1 g{z)) T we obtain 
the following coupled first-order differential equations: 



i d 
a dz 



n 



2 



f{z) 



i d 
a dz 



g(z) = 

' i g(z) = o. 



where Xj = abj/2. The corresponding solution for g(z) 
is easily obtained. We thus obtain the full solution for 
ip(z) shown in Eqs. (Ill) and (pi). 



APPENDIX B: SOLUTION OF THE BOUNDARY 
CONDITIONS TO OBTAIN AN EXCITATION 
SPECTRUM 



From these coupled equations we can derive a second or- 
der differential equation for /(z): 



/"(*) + a 2 



$2(2) 



= ^(f(z)-g(z))S((z-n) + \) 
+ ^{f{z)+g{z))5{{z-n)-\) 



m 

(A4) 



with delta functions at the boundaries and at the center 



of the unit cell. In the regions 



< z 



n < j and 



< z — n < 4, the second order differential equation for 



m is 

/"(z)+a 2 
= 0. 



m 

(A5) 



Since this differential equation is periodic in z, we can 
solve it within a unit cell and use Bloch's theorem to 
extend the solution over all of z. Since < &i(z) has a pe- 
riod of 1/2 and $2(2) has a period of 1, we divide the 
unit cell into two regions (taking n — for simplicity): 



■ ~ < z < i and \ < z < s 



and 



$x(z) 



$2(2) = 



4. In these two regions 

2z, -1/4 < z < 1/4 
2z-l, l/4<z<3/4 



+ 1/2, -1/4<z<1/4 
-1/2, l/4<z<3/4 



For convenience we rewrite Eqs. (|l7|) and (|l^) as 



to (z)- ( A n f+{z) + B n fr{z) 
W) \-A n g+{z) + B n g^{z) 



(Bl) 



for n — j < z < n + j and 



^ 1 ' \-C n g+{z)+D n g 2 [z) ) ' ^ ; 

for n + i < z < n + |. 

Acceptable solutions must be continuous at the inte- 
rior point, z = n + and at the boundaries of the unit 
cell, z = n — j. At z = n + 7 the boundary condition is 

A A + (n + |)+B„ /r(n+i) 
-A, 3+(n+i)+B n 5l (n+i) 

_ f C n f+{n + \) + D n f-{n+\) 



-C n g+(n+\)+D n g 2 (n+±) 
which can be rewritten as 



(B3) 



At z = n — -1 the boundary condition is 



lt{n-\) A-(n-|) \ /A„ 
-ff^n-l) 3x-(n-|) ) \B n 



/ 2 + (n + |) /a"(»+!) \ fC n _i 



-9t(ri+l) g 2 (n+|) / V 



or 



Taking region las — \< z < \ and region 2 as i < z < 
we can write 



(B4) 







where we have used the periodicity of <£i{z) on the right 
hand side of the above equation. From Eqs. (B3) and 
(B4) we can write 



Cn-l 



where 



M 3 " 1 Mi M 3 _1 M 4 . 



(B5) 



(B6) 



Since the Hamiltonian is periodic the eigenvalues p of P, 
an operator which induces a translation of one period, 
must satisfy the Bloch condition: 



(B7) 



where — A < q y < A . The eigenvalues of P are the roots 
of the characteristic equation 



p 2 ~ptr{P) + |P| = 0. 

Clearly 

P± = \ (tr(P) ± VMP)) 2 - 4 
which implies that 
P++P- 



Since q y is real 



cos(2irq y ) = - fr(P). 



Im{ir(P)} = 



(B8) 



(B9) 



(BIO) 



and 



Re{£r(P)} = 2cos(27rg y ). 



(B12) 



The energy spectrum can now be directly calculated us- 
ing this expression. 
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FIG. 1. The square-vortex lattice, showing the A and B 
sublattices and the corresponding unit cell. The edges of the 
unit cell are aligned with the x and y axes that are the nodal 
directions. 
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FIG. 4. A comparison of the 1-D plane-wave expansion (61 
RLVs) (o) with the 2-D plane- wave expansion (33x33 RLVs) 
results (o) for oto — 14 




FIG. 2. Constant-energy contours of Ho and MBZ bound- 
aries of the A and B square sublattices at an — 5. 
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FIG. 5. As in Fig. § for a D = 20 
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FIG. 3. A comparison of the 1-D analytical spectrum along 
the k x axis (o) with the numerical 1-D plane-wave expansion 
(61 RLVs) results (■) for ato = 7. The numerical 1-D results 
are also shown for the spectrum along the k y axis. 
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FIG. 6. The energy spectrum for an — 14 and N y = 41 
and N x = 1(A), 5(x), 9(«), 13(*), 21(»), and 29(o). 
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FIG. 10. The zero bias TDOS for the 1-D plane-wave ex- 
pansion at Qfo = 20 
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^ FIG. 8. The contribution to the TDOS from the nodes at 
k = (±&f,0) at three different energies for an = 14. The 
TDOS is normalized as in Eq. dla). d/2 is the separation, in 
the y direction, between lines ofvortices. Note the shoulders 
forming on either side of the peaks as the energy is increased. 
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FIG. 11. A fit of the magnetic field dependence of the vor- 
tex-core radius as determined from muon spin resonance to 
the l/-\/B (Ref. 25) scaling expected from our analysis. 
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FIG. 9. As in Fig. | for = 20. Note the now very dis- 
tinct shoulders which nave formed on either side of the peaks. 
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FIG. 12. The total DOS in units of N (Tiv F /£) for the node 
k — (kp,0) and an = 14, scaled to show the correspondence 
with the SC calculation for the square-vortex lattice (solid 
line) and in the single- vortex approximation (dashed line). 
Note that both axes scale as yB. The energy is in units of 
E v = Tivf/1- Also shown (thick solid line) is the "averaged" 
quantum DOS, broadened with a Gaussian of full width 0.08 
E v . The inset shows the low energy DOS compared to the SC 
approximations. The averaged quantum DOS is not shown in 
the inset. 
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FIG. 14. The specific heat C v /T for a D = 11 (solid line), 
an = 14 (dotted line), q_d = 17 (dashed line), and an = 20 
(dash-dot line). The inset shows a magnification near T — 
of the same. 
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FIG. 13. As in Fig. |l| for a D = 20. 
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FIG. 15. The specific heat, scaled as 1/ao to show the cor- 
respondence with Cv/T calculated from the Doppler-shifted 
energy spectrum. The inset shows the specific heat with the 
zero magnetic field value subtracted as is done in Ref. 26. 
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